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FOREWORD 


With  the  increasing  use  of  composites  in  aero-space  structure 
components,  survivability  and  vulnerability  of  the  composites  when 
exposed  to  intense  thermal  irradiation  become  a  critical  concern.  New 
materials  with  desirable  surface  characteristics  are  being  explored  to 
mitigate  the  likely  encountered  thermal  loadings;  and  parallel  to  these 
developments,  more  reliable  methods  of  analyzing  the  effects  of  inspired 
surface  heat  flux  on  the  interior  temperature  responses  of  composite 
structures  are  needed. 

The  investigation  reported  here  is  a  preliminary  step  to  examine 
the  basic  needs  for  analytic  tools  to  assess  the  temperature  dis¬ 
tributions  in  composites.  Contained  in  this' report  are  a  new  method  of 
analyzing  cyclic  surface  temperature  fluctuations,  temperature  response 
results  of  composite  panels  due  to  spot-heating,  and  a  comparative  study 
of  the  effect  of  surface  coating  on  the  temperature  rise  of  a  substrate 
layer . 

The  work  was  sponsored  by  Air  Force  Flight  Dynamic  Laboratory  with 
Mr.  Nelson  Wolf  as  the  Technical  Monitor;  it  was  begun  in  June  1983  and 
completed  in  May  1984. 
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1.  INTRODUCTION 


Light  weight,  high  strength  and  dimensional  stability  are  among  the 
principal  advantages  of  composite  materials  in  comparison  to  metals,  and 
consequently  have  led  to  their  ever-increasing  use  in  aerospace  struc¬ 
tures,  For  the  same  reasons,  a  lengthening  list  of  new  composite 
materials  is  emerging  into  the  marketplace  for  specific  applications. 
In  the  new  applications  there  exist  environmental  factors,  heretofore 
uncontemplated,  where  thermal  considerations  are  a  prime  concern. 
Severity  of  the  thermal  environments,  in  which  composite  materials  are 
expected  to  function,  plays  a  critical  role  not  only  in  their  selection, 
but  also  in  the  conceptual  stage  of  composites’  development  as  aerospace 
structural  components.  In  order  to  assess  the  ability  of  composite 
materials  from  a  thermal  viewpoint,  an  analytical  determination  of  the 
temperature  responses  of  candidate  structures  is  of  course  a  first  step 
leading  to,  and  pinpointing,  further  refinements  and  subsequent  develop¬ 
mental  efforts. 

For  thermal  analysis,  the  most  fundamental  characteristic  of 
composite  materials  is  the  effective  thermal  conductivity;  it  is  essen¬ 
tially  reflective  of  the  physical  relationships  between  different  phases 
of  materials  in  a  composite  medium.  In  the  case  of  a  typical  composite 
—  graphite  fibers  in  an  epoxy  binding  matrix  —  the  dispersion  pattern, 
relative  size,  and  density  of  fibers  are  the  governing  criteria  that 
determine  an  effective  conduction  coefficient.  Experimental  effective 
thermal  conductivities  exist  only  for  a  very  few  composite  materials. 
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and  moreover,  they  are  mostly  for  fabricated  specimens  of  fixed  compo¬ 
sitions.  These  ad  hoc  data,  scattered  and  scanty,  are  difficult  to 
organize  and  to  relate  to  one  another  in  order  to  form  engineering 
correlations;  the  situation  is  not  unlike  that  of  the  newly  developed 
aluminum  alloys  decades  ago  when  industry-wide  standards  were  still  in 
their  formative  stages. 

Coupled  with  the  need  for  more  reliable  thermo-physical  properties, 
consideration  must  be  given  to  the  types  of  analysis  for  examining 
thermal  responses  of  composite  materials  in  high-temperature  environ¬ 
ments.  Instead  of  an  overview  of  the  entire  spectrum  of  the  heat 
conduction  phenomena,  a  narrower  but  pragmatic  perspective  of  surviva¬ 
bility  and  vulnerability  of  composite  materials  in  aerospace  structure 
is  adopted.  From  such  a  perspective,  this  report  addresses  two  major 
problem  areas  which  constitute  two  main  critical  tests  for  a  composite 
material  to  survive;  and  they  are  undertaken  in  this  report  not  as 
exhaustive  accounts  from  an  operational  viewpoint,  but  as  demonstrative 
studies  for  establishing  the  methodology  of  each.  Results  of  the  two 
case  studies  are  discussed  in  terms  of  thermo-physical  properties  of  the 
composites  and  their  parametric  relations  with  calculated  temperature 
responses.  From  the  calculated  results,  guidelines  are  established  for 
an  assessment  of  composites1  ability  to  withstand  the  intended  thermal 
environment . 

The  first  problem  concerns  the  localized  thermal  heating  of  a  panel 
by  a  cylindrical  beam  of  irradiation.  As  composite  materials  in  general 
are  characterized  by  orthotropical  conductivities,  the  analysis  focuses 
on  the  effect  of  the  in-plane  thermal  conductivity  upon  the  heat 
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penetration  pattern.  The  second  problem  deals  with  the  temperature 
response  of  a  multi-layer  composite  subjected  to  cyclic  surface  heating, 
which  is  a  representative  encounter  for  an  orbiting  space  vehicle.  Both 
problems  are  limited  to  their  parametric  performances  of  temperature 
responses,  and  their  subsequent  physical  phenomena  —  melting  or  abla¬ 
tion  in  the  first  case  and  buckling  with  alternating  thermal  stresses  in 
the  second  —  are  not  discussed  in  this  report. 


2.  AXIS-SYMMETRICAL  SPOT-HEATING  OF  COMPOSITE  SLABS 


Consider  a  slab  with  a  thermal  irradiation  loading  by  a  cylindrical 
beam  of  constant  heat  flux.  The  material  of  the  slab  is  such  that  the 
thermal  conductivity  in  the  plane  of  the  slab  is  isotropic  but  differs 
from  the  conductivity  in  the  depth  direction.  Such  a  combination  of 
thermal  conductivities  typifies  composites  with  fibers  oriented  in  the 
plane  of  the  slab.  Composites  with  fibers  in  overlay  patterns  of  0/90, 
0/±45/90  all  fall  into  this  category. 

The  slab  is  initially  at  a  uniform  temperature,  arbitrarily  taken 
to  be  zero,  and  as  spot-heating  proceeds,  the  temperature  rise  of  the 
panel  is  to  be  determined. 

2.1  Analysis .  Figure  1  depicts  schematically  the  thermal  system 
under  consideration.  Thermal  irradiation  is  confined  to  a  radius  of  a 
and  is  of  intensity  Q.  Axis-symmetry  is  assumed  and  the  heat  diffusion 
equation  together  with  its  boundary  conditions  are: 
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Figure  1.  Cylindrical  Irradiation  Beam  Configuration 

-V - 
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Z  Figure  2.  Node  Definition 
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at  z  «  0,  r  <  a 


(2) 


z  -  0,  r  >a;  and  at  z  =  w 


(3) 


Except  at  the  irradiation  beam,  all  surfaces  are  insulated  and  the  slab 
is  considered  infinite  in  the  in-plane  direction.  The  characteristic 
parameters  are  a,  w,  Q  and  the  thermal  properties  (PC),  k^,  and  k^,  all 
of  which  are  needed  to  express  the  results  in  non-dimensional  forms. 

The  terms  of  the  right-hand  side  of  equation  (1)  represent  of  course  the 
heat  fluxes  in  the  r  and  z-direction  respectively,  which  lead  to  the 
following  non-dimensional  variables: 


(4) 


z  =  z/a 


(5) 


e  =  (a  e/a2) 

Z 


(6) 


where  a is  the  thermal  diffusivity  and  is  given  by 


a  =  k  /(pc) 


(7) 


z  z 


The  preceding  transformations,  particularly  equations  (4)  and  (5), 
indicate  a  scale  change  in  the  r-direction  because  of  different  thermal 
conductivities  k^  and  k^.  By  defining  a  non-dimensional  temperature  as 
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T  - 


(8) 


the  preceding  equations  are  rendered  non-dimensional  as  follows: 


9T  =  82T  JL  3T  82T 

90  9r  r  9r  9z 


(9) 


The  boundary  conditions  become 


9T 

9z 


-1 


z 


0, 


(10) 


—  =  0  z  =  0,r>/kz/k^  (11) 

9z  ^ 


-Z  -  0  z  =  w/a 

9z 


(12) 


Initially,  before  heat  flux  at  z  =  0,  the  temperature  is  zero,  or 


T  =  0  0  =  0 


(13) 


The  different  scale  changes  in  the  r  and  z  directions  are  made  necessary 
because  of  different  thermal  conductivities  in  these  two  directions. 
Solution  of  the  preceding  set  of  equations  can,  in  principle,  be  ob¬ 
tained  by  a  rigorous,  exact  approach,  for  example,  by  a  Laplace  trans¬ 
form  in  0,  a  Fourier  transform  in  z  and  a  Hankel  transform  in  r.  Such 
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an  approach  is  conceptually  satisfying;  however,  the  inversion  proce¬ 
dures  necessary  for  computation  may  become  so  involved  that  the  method 
is  itself  rendered  impractical  and  moot.  Hence,  in  order  to  produce 
numerical  results  for  making  engineering  judgements,  the  finite- 
difference  approach  is  used  and  its  implementation  is  described  next. 

2.2  The  Finite-Difference  Solution.  For  diffusion  problems,  the 
simplest  scheme  is  a  combination  of  three-point  central  difference  in 
space  and  one  time-step  marching  using  explicit  algebra.  It  however 
suffers  from  the  disadvantage  of  being  unstable  which  can  only  be  cured 
by  using  smaller  time  steps.  Implicit  schemes  are  known  to  be  uncondi¬ 
tionally  stable  but  requires  the  inversion  of  large  matrixes.  In 
one-dimensional  problems,  the  matrixes  are  tri-diagonal;  but  in  multi¬ 
dimensional  ones,  the  matrixes  are  not  so  simple.  And  so  far,  no 
simple  algorithms  have  been  devised. 

Instead  of  the  fully  implicit  and  fully  explicit  method,  a  viable 
choice  is  the  alternating-direction  explicit  method  [1]  which  has  a 
greater  degree  of  stability  and  is  therefore  less  time  consuming  than 
the  explicit  method.  In  addition,  the  necessary  algebraic  manipulation 
is  much  less  involved  than  the  implicit  method,  especially  for  multi¬ 
dimensional  problems . 

Since  the  heat  beam  radius  is  a  governing  dimension  of  the  problem, 
the  size  of  a  finite-difference  grid  begins  at  the  heat  spot.  Along  its 
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radius,  5  divisions  are  used.  The  exact  numerical  value  of  Ar 

is  dependent,  of  course,  on  the  relative  magnitudes  of  and  k^ .  By 

definition  ?  =  (r/a)  /k  /k  ,  hence  from  r  =  0  to  r  =  a  there  are  five 

\j  z  r 

equal  increments,  each  equal  to  Ar  =  0.2/k  /k  ,  with  the  ratio  of  the 
two  conductivity  values  varying  from  one  problem  specification  to 
another.  To  attain  equal  accuracy  in  the  z-direction,  equation  (9) 
requires  that  Az  =Ar.  In  this  way  the  mesh  size  is  made  flexible, 
depending  on  the  relative  magnitudes  of  the  conductivities.  The  total 
number  of  nodes  in  the  zf-direction  is  determined  by  the  slab  thickness 
(w/a) . 

It  should  be  mentioned  that  in  the  finite-difference  algorithm  usecb 
for  the  two-dimensional  problems,  a  single-indexed  temperature  array  is 
used.  Customarily,  a  double-indexed  array  is  for  keeping  track  of  the 
nodal  positions.  However,  owing  to  the  undetermined  numbers  of  modes  in 
both  directions,  single-indexed  arrays  are  better  suited  to  the  situa¬ 
tion;  the  node  positions  can  be  easily  identified  by  the  array  index  in 
question.  For  example,  if  there  are  20  nodes  in  ~z  and  500  in  ~r9  a 
single  array  of  variables  with  index  from  1  to  10000  is  needed.  The 
same  array  can  also  be  used  for  500  nodes  in  z"  and  20  in  r.  A  simple 
accounting  of  the  index  value  serves  to  identify  the  node  in  question. 

In  presenting  the  various  finite-difference  formulas,  however,  the 
usual  two-index  notations  are  used  for  simplicity  and  ease  of  dis¬ 
cussion. 


*  The  effects  of  using  10  divisions  are  documented  in  Appendix  A 
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The  Finite-Difference  Formula  for  Regular  Nodes.  A  node  located  in 


the  interior  of  the  slab  but  not  the  axis  of  symmetry  is  termed  regular 

node.  Using  two-indexed  notations,  T  is  the  temperature  at  the  node 

i » J 

position  identified  by  the  two  indexes.  On  the  exposed  surface  where 
heating  occurs,  i  =  1  and  on  the  axis  of  symmetry,  j  “  1;  nodes  on  these 
locations  are  not  regular  nodes.  The  finite-difference  formula  for 
regular  nodes  is  based  on  a  r-z  mesh  of  equal  increments  as  are  for 
other  nodes.  A  five-point  cluster  of  nodes  is  illustrated  in  Figure  2 
together  with  their  definitions. 

According  to  the  ADE  method,  the  second-order  derivative  on  the 
right  of  equation  (9)  is  replaced  by  a  time-splitting  procedure.  Thus 


32T 


[(ii, 


j+i 


-  t  .  . )  - 
iij 


(Tn  -  Tn  )1  h7 

Ui,j  ‘i.J-l/J/ 


(14) 


Superscript  n  denotes  the  nodal  values  at  a  later  time  step,  or  "new" 
values;  otherwise,  current  values  are  meant.  The  same  increment  of  Ar 
and  Az  is  denoted  by  a  common  symbol  A.  Similarly,  the  following  is 
also  valid: 
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9z 


(15) 


The  first-order  derivative  of  equation  (9)  is  given  by  a  central- 
difference  approximation,  again  using  split-time  values. 
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(Ti,j+l  -  tL-i)/(24) 


(16) 


3T 

dT 


3T 

30 


(17) 


where  6  represents  the  time  step  6=A9.  Incorporating  these  individual 

expressions  into  equation  (9),  the  resulting  formula  for  the  calculation 

of  Tn  .  is  given  by 
1 » J 


+  <«M2)  [TtjJ+1  -  + 

/[l  +  2  (6/A2)] 


(18) 


Elementary  considerations  of  the  stability  question  lead  to  a  positive 

coefficient  of  T  .  on  the  right  of  equation  (18).  Thus  the  time  step  6 
1  >  J 

is  related  to  the  space  step  A  by 


6  <  ( A 2 /2 )  09) 

In  fact,  as  shall  be  developed  later,  the  time  step  6  is  taken  to  be 

6  =  A2 /2  (?°) 

Even  though  the  right-hand  side  of  equation  (18)  contains  nodal  values 
at  a  new  time,  these  are  already  known  from  the  preceding  calculations 
if  the  finite-difference  computations  are  proceeding  in  the  direction  of 
increasing  index  numbers. 
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The  Finite-Dif ference  Formula  for  Central  Nodes  (r  «  0) .  For  the 
nodes  located  on  the  axis  of  symmetry,  where  r  »  0,  equation  (18)  is  not 
valid  since  it  contains  a  factor  of  1/r  on  its  right-hand  side.  The 
axis-symmetrical  nature  of  the  temperature  distribution  can  be  repre¬ 
sented  by  an  algebraic  expression. 


T  =  Tj  +  br2 


(21) 


which  satisfies  the  requirement  3T/3r  c  0  at  r  =  0,  The  temperature  at 


r  =  0  is  denoted  by  T^;  at  r  *  A ,  the  temperature  denoted  by  is 


related  to  the  coefficient  b  in  equation  (21).  Equation  (21)  then 
becomes 


T  =  Tj  +  (T2  -  T1)(r2/A-2) 


(22) 


which  is  valid  for  small  values  of  r.  The  r-derivatives  of  equation  (9) 
can  then  be  obtained  from  equation  (22)  and  they  become 


(23) 


(24) 


Expressed  in  terras  of  two-indexed  notations  and  using  the  split-time 
schedule,  the  sum  of  these  two  terms  is  given  by 


(25) 
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For  the  z-direction,  the  split-time  scheme  results  in 


ill 

dz2 


[Ti+i,i 


(26) 


From  equations  (25)  and  (26),  the  finite-difference  formula  for  central 
nodes  is 

TM  ’  [If/S’I|T1«,1  +  Ti-I.i  +  4V'  +  V11  -  «/A2)l] 

/[  1  +  5 (6/A2) ]  (27) 


Stability  considerations  of  equation  (27)  in  calculating  the  new  temper¬ 
ature  lead  to  the  criterion 


6  <  A2 


(28) 


which  makes  the  coefficient  of  T  .  of  equation  (27)  positive. 

i ,  I 

The  criteria  of  equations  (19)  and  (28)  are  both  required  and 
obviously  the  former  is  more  stringent;  hence  the  time  step  6  is 
governed  by  that  of  equation  (20). 

The  Finite-Difference  Formula  for  Boundary  Nodes.  The  boundary 
condition  at  the  slab  top,  "z  *=  0,  is  described  through  the  first  deriva¬ 
tives  of  the  temperature  by  equation  (10).  Starting  from  the  top 
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surface,  the  z-index  is  i  =  1,  2,  etc.  The  first  derivative  by  a 
3-point  formula  is  therefore 


*  (4  T  -  T,  .  -  3T  .)/ 2A  (29) 

^  > 3  3 ,  j  1,  j 


For  nodes  located  in  the  heating  beam,  r<a,  equation  (10)  governs  which 

raises  the  surface  temperature  to  a  new  high  value  in  successive 

time  steps.  By  setting  the  derivative  of  equation  (29)  at  -1  and 

replacing  T  .  by  T*}  ,  a  finite-difference  formula  for  the  new  tempera- 

1  >  J  1  y  J 

ture  is  then  established: 


(4  T2,J  -  T3,j  +  2A)/3 


(30) 


For  nodes  lying  outside  the  irradiation  beam,  the  condition  of  zero  heat 
transfer  results  in 


l.j 


(4  T 


2>i 


-  T3,j)/3 


(31) 


2.3  Scope  of  Computations 

As  this  investigation  is  intended  to  exemplify  and  demonstrate  the 
need  of  major  areas  of  further  research  and  to  establish  principal 
thermal  criteria  for  composite  application  in  aero-space  structure, 
computational  efforts  were  limited  to  the  variations  of  a  few  major 
parameters  so  that  conclusions  can  be  drawn  to  pinpoint  the  future 
development  of  composite  materials  from  a  thermal  protection  viewpoint. 
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Accordingly,  starting  with  the  reference  case  of  a  panel  with  isotropic 
conductivity,  i.e.  ,  two  additional  conductivity  variations  were 
analyzed:  they  are  (i)  radial  (in-plane)  thermal  conductivity  four 
times  as  large  as  the  transverse  value,  k^  =  4k^ ,  and  (ii)  k^  *  9k^  as 
the  extreme  case.  Altogether,  the  conductivity  ratios  are  1,  4,  and  9; 
these  are  considered  adequate  or  at  least  representative  of  the  expected 
variation . 

Another  geometrical  parameter  is  the  ratio  of  the  panel  thickness  w 
to  the  spot  radius  a.  In  the  computational  effort  of  this  analysis, 
values  of  (w/a)  of  1,  2,  and  4  were  taken.  Hence  altogether,  there  are 
nine  combinations  of  conductivity  ratios  and  panel  thickness  ratios. 
From  these  combinations  of  the  parameters,  computed  temperature 
responses  due  to  impulsive  spot  heating  are  examined  to  formulate 
performance  criteria  for  composite  materials  serving  as  thermal  barriers 
to  surface  irradiation. 

From  a  survivability  and  vulnerability  standpoint,  the  criteria  are 
temperature  rise  of  the  heating  spot  as  heating  proceeds,  depth  of  heat 
penetration  into  the  substrate,  and  back  surface  temperature  rise.  The 
temperature  rises  in  conjunction  with  other  relevant  considerations  — 
such  as  glass /transition  temperature  or  melting  temperature,  localized 
buckling  and  de-lamination,  —  form  the  basis  of  selection  and  applica¬ 
tion  of  composite  materials. 

2.4  Results  and  Discussion 

Spot-Surface  Temperature  Rise.  The  temperature  rise  at  the  center 
of  the  irradiation  beam  is  the  most  important  parameter  of  the  thermal 
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response  for  judging  how  well  a  composite  material  endures  a  severe 
environment.  In  general  terms,  high  temperature  as  a  result  of  external 
thermal  loads  engenders  various  damages  which  may  lead  to  progress 
incapacitation  of  the  aerospace  structures:  possible  damages  include 
localized  melting;  structural  buckling  and  fracture;  layer  delamination 
and  crack  initiative  or  enlargement;  and,  of  course,  thermal  puncture 
due  to  mass  removal.  Multitude  of  these  damage  mechanisms  precludes  a 
complete  analysis  of  each;  instead,  emphasis  is  placed  on  the  role  of 
the  in-plane  thermal  conductivity  in  affecting  the  temperature  response 
due  to  a  confined  irradiation  beam. 

Starting  with  the  reference  case  of  isotropic  conductivity,  i.e., 
k^_  =  anc*  ^°r  Pane^s  thickness  ratios  of  w/a  =1,2,  and  4,  the 

spot  center  temperature  rises  are  shown  in  Figure  3.  The  back  surface 
temperature  rises  are  also  indicated.  The  frontal  spot  surface  tempera¬ 
tures  are  identified  by  "s"  and  the  back  center  spot  temperatures  by  "b" 
in  Figure  3  with  the  numerals  preceding  s  or  b  indicating  the  slab 
thickness  ratio. 

An  obvious  trend  of  the  temperature  curves  is  that  the  magnitudes 
are  mitigated  by  increased  thickness  from  curve  nos.  Is,  2s,  to  4s; 
temperature  reduction  is  more  pronounced  for  the  back  surface  tempera¬ 
ture.  These  numerical  values  establish  a  reference  with  which  other 
calculations  can  be  compared.  The  reference  performances  in  Figure  3 
also  illustrate  the  thermal  protection  available  by  straight¬ 
forward  thickening  of  thermal  shields;  in  this  case  the  additional 
panelty  to  reduce  the  front  and  back  surface  temperature  rise  is  in 
direct  proportion  to  the  added  thickness. 
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Figure  3.  Heat  Spot  and  Back  Temperatures  for  Isotropic 
Conductivity,  k^/k^  =  1 


-16- 


-  ( kz  To  /qQ  ) 


1.2 


.0  - 


l  =  (kr/k2  ) 


0.8 


0.6 


W-9 


0.4 


0.2- 


0I _ I _ I _ I _ I _ L 

0  2  4  6  8  10 

e=(aze/Qz) 

Figure  4.  Heat  Spot  Temperature  Rises  for  Large  In-Plane 
Thermal  Conductivities,  w/a  =  1 
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Figure  5.  Heat  Spot  and  Back  Temperature  Rises  by  Thermal 
Shieldings,  k^/k^  =  4 
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Figure  6.  Heat  Spot  and  Back  Temperature  Rises  by  Thermal 
Shieldings,  kr/kz  =  9 
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It  is  significant  to  note  that  as  heating  proceeds,  the  front 
surface  and  the  back  surface  temperatures  tend  toward  a  fixed  differen¬ 
tial.  Pairs  of  curves,  for  example.  Is,  and  lb,  gradually  become 
parallel  to  each  other.  This  phenomenon  indicates  that  within  the  slab, 
the  net  heat  flow  occurs  in  the  radial  direction  and  is  essentially 
spreading  along  the  slab  plane.  In  the  early  part  of  the  heating 
period,  say  0<1,  the  heat  spot  temperature  rises  are  similar  to  those 
of  an  infinite  half-space  with  an  impulsive  surface  heat  flux.  As 
heating  proceeds,  radial  heat  spread  becomes  effective  and  ultimately 
dominant  in  limiting  the  temperature  change,  as  evidenced  by  the  bending 
of  the  temperature  curves. 

This  observation  of  the  temperature  trends  lends  itself  to  the 
interpretation,  and  indeed,  expectation,  that  composite  panels  with 
large  in-plane  thermal  conductivities  would  promote  rapid  lateral 
temperature  spread  and  thereby  reduce  the  front  and  back  surface  temper¬ 
ature  rises. 

To  assess  quantitatively  the  effect  of  large  in-plane  conductivity, 
additional  computations  are  undertaken. 

For  a  slab  with  a  thickness  ratio  of  w/a  **  1  and  conductivity 
ratios  of  k^/k^  =  1»  4,  and  9,  the  computed  temperature  responses  are 
grouped  together  in  Figure  4.  It  is  readily  apparent  that  the  tempera¬ 
ture  curves  for  k^/k^  *  4  and  9  show  very  early  trends  toward  the 
phenomenon  of  asymptotic  radial  heat  flow  from  the  heat  source  of  the 
irradiation  beam.  The  surface  spot  temperature  rises  are  thus  very  much 
reduced  compared  to  the  case  of  k^/k^  “  moreover,  the  temperature-rise 
limiting  feature  by  virtue  of  the  higher  in-plane  conductivity  is  more 
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effective  than  simple  thickening  of  the  thermal  shield  of  an  isotropic 
material,  as  shown  in  Figure  3. 

To  further  bring  out  the  relative  performances.  Figure  5  groups 
together  the  temperature  calculations  for  k^/k^  =  4  and  slabs  of  thick¬ 
ness  ratio  of  1,  2  and  4.  Figure  5  is  a  parallel  to  Figure  3;  the 
latter  is  for  k^/k^  =  1*  Similarly,  Figure  6  contains  the  temperature 

data  for  k  /k  =9  and  the  three  thickness  ratios.  The  performance  data 
r  z 

in  these  two  figures  demonstrate  that  large  in-plane  conductivities  lead 

to  an  early  (in  time)  asymptotic  heat  flow  pattern  of  radial  spread, 

thus  reducing  the  front  and  back  surface  temperature  rises.  Further 

reductions  of  the  temperatures  are  possible  with  additional  thermal 

shield  thickening,  but  the  return  is  diminishing  for  higher  ratios  of 

k  /k  . 
r  z 

Asymptotic  Trend  of  Temperature  Variations.  As  shown  by  the 
preceding  discussion,  the  heat  spot  temperature  increases  at  first 
according  to  that  of  a  semi-infinite  medium  with  an  impulsive  heat  flux 
across  the  entire  surface.  After  a  short  time  interval,  the  finite 
width  of  the  slab  and  the  limited  heat  beam  radius  begin  to  affect  the 
temperature  rise  pattern:  The  pattern  eventually  becomes  heat  con¬ 

duction  In  a  finite-thickness  sheet  with  heat  input  on  the  internal 
boundary  of  a  cylindrical  surface.  Such  an  asymptotic  configuration  is 
defined  in  Figure  7,  showing  an  internal  boundary  of  radius  a  equal  to 
the  heat-beam  radius  and  a  heat  flux  magnitude  q^. 

Such  a  heat  flow  problem  has  been  analyzed  and  presented  in  [2]. 
For  large  times,  the  temperature  rise  at  the  heating  boundary  r  =  a  is 
given  by 
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2k  T 
r  a 


(32) 


=  Log  [4(0  0/a2)]  -  0.577 
aqr  6e  r 

A  significant  conclusion  of  equation  (32)  is  that  at  large  times, 
the  temperature  rise  varies  according  to  the  log  of  time  0  .  To  test 
such  a  correlation,  the  temperature  rises  computed  for  all  nine  cases 
are  grouped  together  and  plotted  vs.  the  log  of  the  lapse  time  0  . 

Since  the  principal  variable  in  the  case  studies  is  the  in-plane 
thermal  conductivity  k^,  and  the  non-dimensional  temperature  and  time 
use  and  a as  reference  quantities  in  the  ordinate  and  abscissa  of 
Figure  8,  in  this  way  the  relative  magnitudes  of  T  give  a  direct  in¬ 
dication  of  the  surface  heat  spot  temperatures. 

To  aid  identification,  these  nine  performance  curves  are  numbered. 
The  first  set,  nos.  1,  2,  and  3,  shows  the  temperature-time  histories 

for  the  three  slab  thicknesses  w/a  ~  1,  2  and  4,  and  k  /k  =  1 ;  the  next 

r  z 

set,  nos.  4,  5,  and  6,  is  for  k  /k  =*4  and  the  same  thickness  varia- 

r  z 

tions;  and  the  last  set,  nos,  7,  8,  and  9,  for  k^/k^  «  9.  Having  all 
nine  cases  together  in  a  single  graph  again  demonstrates  the  alleviation 
of  the  heat  spot  temperature  rise  by  virtue  of  higher  in-plane  thermal 
conductivities,  as  contrasted  to  simple  thicker  thermal  shields.  More 
important  is  the  feature  that  for  longer  heating  time  0  >  4,  the  semi- 
logrithmic  plot  in  Figure  8  displays  linear  relations,  thus  essentially 
confirming  the  trend  implicit  in  equation  (32).  The  fact  that  for  each 
case  shown  in  Figure  8,  the  slope  of  its  asymptotic  variation  is 
different  from  the  others  can  be  reconciled  by  noting  the  left-hand  side 
of  equation  (32)  and  the  non-dimensional  ordinate  of  Figure  8.  Equation 
(32)  contains  a  radial  heat  flux  q^,  as  is  indicated  in  Figure  7, 
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which  is  related  to  the  surface  heat  flux  Q  in  the  irradiation  beam  by 
the  following  equivalence: 


Tra2Q  =  2  Trawq^ 

If  q^  is  eliminated  from  equation  (32)  by  using  the  above  relation, 
equation  (32)  is  converted  to 

?  2 

(wk  T  /a  Q)  =  0.25  Log  (a  6/a  )  +  constant  (33) 

r  a  ^  °e  z 

where  denotes  the  temperature  of  the  slab  at  the  beam  edge  r  =  a. 

The  re-organized  temperature  parameter  on  the  left-hand  side  of  equation 
(33)  is  then  used  to  correlate  the  calculated  temperature  rise  at  the 
heat  spot  center  —  but  not  the  beam  edge.  The  resulting  relations  for 
the  nine  cases  are  shown  in  Figure  9,  where  a  prominent  feature  is  that 
the  asymptotic  slopes  for  these  curves  are  nearly  equal  to  0.30,  which 
results  in  an  equation  of  the  form, 

2  2 

(wk  T  /a  Q)  =  0.30  Log  (a  0/a  )  +  constant  (34) 

r  o  e  z  v  ' 

where  T^  is  the  heat  spot  center  temperature. 

Equations  (33)  and  (34)  differ  from  each  other  in  the  coefficients 
of  0.25  and  0.30,  which  are  attributable  to  the  different  locations  in 
the  asymptotic  equivalent  configuration  of  Figure  7,  where  heating 
occurs  at  r  =  a  and  in  Figure  9,  where  T^  at  r  =  0  is  used.  Examination 
of  the  temperature  profiles  for  the  nine  case  studies  indicates  that 
there  is  a  decline  of  the  distribution  from  r  =  0  to  r  =  a  with  a  ratio 
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of  T  /T  =  1.25;  in  other  words,  if  T  is  used  in  lieu  of  T  in  Figure 
o  a  a  o 

9,  the  resulting  slope  would  become  0.30/1.25  =  0.24,  which  is  in 
excellent  agreement  with  the  asymptotic  radial  heat  spread  requirement. 

It  should  be  emphasized  that  the  matching  of  the  asymptotic  slopes 
with  the  simplified  configuration  of  Figure  7  is  not  complete,  for 
Figure  9  indicates  that  these  variations  still  differ  from  one  another 
by  a  constant  term  in  equation  (34);  i.e.,  the  curves  are  at  different 
levels.  This  must  be  attributed  to  the  different  combinations  of 
(k^/k^)  and  (w/a)  in  the  finite-difference  calculations;  each  one  has  a 
different  temperature  vs.  time  history  in  terms  of  heat  storage,  etc. 
before  reaching  their  respective  asymptotic  states. 

A  factor  of  significance  is  that  if  the  temperature  response  of  a 
composite  panel  is  to  be  analyzed,  calculations  can  be  terminated  at  a 
time  asymptotic  radial  heat  spread  becomes  established,  after  which  time 
a  semi-logarithmic  formula  can  be  used  for  extrapolating  the  spot 
surface  temperature  to  larger  values  for  9  . 

Isotherms  and  Heat  Penetration  Depths.  To  further  delineate  the 
effect  of  the  in-plane  thermal  conductivity,  lines  of  constant  tempera¬ 
ture  extracted  from  the  numerical  results  are  presented.  Since  there 
are  a  number  of  parametric  variations  in  the  cases  analyzed,  isothermal 
contours  at  9=  10  only  are  considered.  The  isotherms  are  identified  in 
terms  of  fractions  of  the  heat  spot  surface  temperature;  this  repre¬ 
sentation  is  necessary,  for  the  individual  temperatures  of  the  heat  spot 
surfaces  are  all  different  from  each  other. 
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Figure  7.  Model  of  Constant  Heating  of  an  Infinite  Medium 
Internally  Bounded  by  a  Cylindrical  Surface 
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Figure  9.  Asymptotic  Variations  of  Heat  Spot  Temperature  Rises 
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Shown  in  Figure  10  are  isothermal  contours  which  define  temperature 

rises  corresponding  to  0.8,  0.6,  and  0.4  of  the  individual  hot  spot 

temperature  rises  for  the  cases  of  k^/k^  =  1 »  4 ,  and  9;  all  for  w/a  =  1. 

For  the  case  of  k  /k  =1  shown  in  Part  (a)  of  Figure  10,  the  tempera- 
r  z 

ture  contours  beyond  the  heating  beam  radius  of  a  are  nearly  transverse 
to  the  slab;  this  is  of  course  indicative  of  near  one-dimensional 
temperature  distributions  in  which  the  distribution  is  asymptotic, 
gradually  approaching  that  of  Figure  7.  For  larger  in-plane  thermal 
conductivities,  k^/k^  83  4  and  9,  temperature  contours  shown  in  Parts  (b) 
and  (c)  of  Figure  10  signify  higher  radial  (in-plane)  spread  than  in  the 
depth  direction.  Not  only  is  the  back  surface  temperature  reduced  by 
rapid  heat  dissipation  near  the  heat-input  surface,  but  more  so  is  the 
temperature  rise  of  the  surface  heat  spot  itself.  At  heating  time  6= 

10,  for  which  Figure  10  is  valid,  the  spot-center  temperatures  for  the 
three  conductivity  ratios  are  1.49,  0.68  and  0.45,  respectively.  This 
observation  is  most  significant  when  the  criterion  of  burn-through  of 
thermal  shields  is  considered. 

Similar  to  the  contours  in  Figure  10,  the  cases  of  w/a  =  4  and 

k  /k  =1  and  9  are  displayed  in  Figures  11  and  12.  For  isotropic 
r  z 

conductivity,  the  contours  (Figure  11)  at  0  =  10  are  very  nearly  circu¬ 
lar  (spherical  in  reality),  as  would  be  expected;  in  contrast,  for  k^/k^ 
=  9,  temperature  rises  are  concentrated  near  the  heat  source  but  spread 
appreciably  along  the  plane,  resulting  in  contours  of  flat  ellipses. 

Another  feature  of  significance  is  the  profile  of  the  slab  tempera¬ 
ture  on  the  heat  input  surface  *z"  *  0;  Figures  13  and  14  show  the  lateral 
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temperature  variations  for  k^/k^  «  1  and  9,  each  with  w/a  =  1  and  4, 
respectively,  in  these  two  illustrations.  As  in  the  case  of  the 
isothermal  contours,  the  temperature  profiles  are  presented  in  terms  of 
their  ratios  to  the  temperature  of  the  heat  spot  center.  In  spite  of 
wide  variations  of  the  parameters  (w/a)  from  1  to  4  and  the  conductivity 
ratio  (k^/k^)  from  1  to  9,  the  surface  temperature  profiles  are  re¬ 
markably  similar  to  each  other:  a  very  nearly  flat  region  within  the 
beam  radius  and  sharp  drop  from  the  beam  edge  outward.  It  is  parti¬ 
cularly  worth  noting  that  the  temperature  at  the  beam  edge  is  approxi¬ 
mately  0.8  of  the  temperature  at  beam  center  and  the  temperature  ratio 
is  reduced  to  0.3  at  a  distance  of  two  radii  from  the  center.  The 
former  relation  is  used  in  correlating  the  surface  temperature  rise  with 
the  asymptotic  heat  flow  configuration  of  Figure  7  and  equation  (34). 

2.3  Summary  and  Recommendations 

At  the  end  of  0  =  10,  the  heat  beam  center  temperatures  and  the 
back  temperatures  for  the  nine  combinations  of  parameters  are  listed  in 
Table  I.  Starting  with  the  case  of  k^/k^  =  1  and  w/a  «  1,  the  heat  beam 
spot  temperature  is  reduced  from  1.49  to  1.015  by  increasing  the  thick¬ 
ness  to  w/a  =  4.  However,  a  greater  reduction  of  the  surface  tempera¬ 
ture  can  be  achieved  by  increasing  the  in-plane  thermal  conductivity: 

Fy  using  k^/k^  =  4  and  9,  the  surface  temperature  becomes  0.682  and 
0.450  respectively.  This  observation  demonstrates  the  need  for  compos¬ 
ites  with  larger  in-plane  thermal  conductivities  than  the  transverse 
values . 
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Aside  from  concluding  that  the  in-plane  thermal  conductivities  of 
composite  labs  have  a  predominant  influence  on  the  thermal  responses  in 
general,  this  investigation  also  demonstrates  other  typical  analyses 
among  many  needed  for  a  survivability  and  vulnerability  evaluation  ol 
protective  materials.  The  types  of  investigations  needed  can  be, 
broadly  speaking,  classified  into  (i)  short-term  damages  characterized 
by  imminent  incapacitation,  such  as  burn-through  melting,  large-deformation 
buckling  and  fracture  and  (ii)  long-term  cumulative  damages,  such  as 
thermally-induced  stress  concentration  around  voids,  cracks,  and  delin¬ 
eation  of  layers,  etc.  Indeed,  the  analysis  of  spot-heating  treating  in 
this  report  —  even  though  limited  to  the  temperature  effect  along  — 
belongs  to  the  first  category,  for  if  the  heating  process  continues, 
burn-through  or  pit-forming  is  the  next  stage  of  occurrence. 
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Figure  11.  Spot  Heating  Isothermal  Contours  for  w/a 
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Figure  12.  Spot  Heating  Isothermal  Contours  for  w/a  =  4  and  k  /k 
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Figure  13,  Surface  Temperature  Profiles  at  Various  Conductivity  Ratios,  w/a 


Front  Surface 


-35- 


Figure  14 .  Surface  Temperature  Profiles  at  Various  Conductivity  Ratios,  w/a 


Table  I.  Calculated  Heat  Spot  and  Back  Temperatures* 


(e=  io) 


<VV 

w/a  =  1 

w/a  =  2 

w/a  =  4 

1 

1.495 

1.121 

1.015 

1.026 

0.395 

0.  107 

4 

0.682 

0.578 

0.550 

0.305 

0.  106 

0.028 

9 

0.450 

0.402 

0.388 

0.  142 

0.048 

0.009 

*  The  top  and  bottom  figures  indicate  the  frontal  and  back  surface 
temperatures  respectively;  the  temperature  is  defined  by 
T  =  k^T/aQ. 
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3.  TRANSIENT  HEAT  FLOW  IN  MULTI-LAYERED  COMPOSITES 


Aside  from  the  axis-symmetrical  spot-heating  problem  analyzed  and 
discussed  in  the  preceding  section,  one-dimensional  transient  heating  of 
multi-layered  composites  constitutes  a  fundamental  class  of  problems  — 
not  necessarily  from  the  viewpoint  of  its  novelty  or  difficulty  —  hut 
from  the  viewpoint  of  providing  realistic  temperature  distributions 
throughout  the  composite  bulk  so  that  the  resultant  thermal  stresses  can 
be  calculated  based  on  more  realistic  thermal  analyses. 

The  first  problem  of  transient  conduction  analyzed  in  this  part  of 
the  report  is  the  heating  by  a  constant  heat  flux  of  a  multi-layered 
composite  panel.  From  a  vulnerability  standpoint,  assessment  of  compos¬ 
ites  used  in  space  structure  must,  by  necessity,  involve  the  constant 
heat  flux  criterion  in  order  to  establish  the  temperature  excursion 
which  a  candidate  composite  structure  may  undergo.  In  this  connection, 
relevant  parameters  are  the  heat  flux  intensity,  its  duration,  and  the 
thermal  properties  of  the  composite  materials.  The  temperature  re¬ 
sponses  derived  therefrom  constitute  the  input  to  the  next  phase  of 
analysis  for  the  resulting  thermal  stresses  and  possible  deformation, 
buckling,  etc.,  which  are  based  on  known  material  thermoelastic  prop¬ 
erties  and  the  geometrical  specification  of  the  space  structural  shape. 

The  second  problem  discussed  in  this  section  pertains  to  the 
periodic  heating  and  cooling  which  comes  about,  for  example,  in  an 
orbiting  space  structure.  During  a  portion  of  its  cycle,  solar  irra¬ 
diation  raises  the  exposed  surface  temperature  and  during  its  stay  in 
the  shadow,  the  surface  temperature  lowers.  The  cyclic  surface 
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temperature  fluctuation  is  a  direct  consequence  of  the  fluctuating 
surface  heat  flux  about  its  cyclic  average  value.  Such  a  periodic  heat 
flux  variation  and  surface  temperature  excursion  are  of  course  mutually 
convertible  and  can  be  expressed  by  a  Fourier  series  whose  basic  period 
is  the  orbiting  time,  with  higher  harmonics  taking  into  account  irreg¬ 
ularities  from  a  sine-wave  representation. 

The  development  herein  focuses  more  on  the  methodology  which  is 
heretofore  unavailable  than  on  the  complex  parameters  involved  in  an 
orbiting  event.  From  an  application  viewpoint,  the  method  of  analysis 
easily  accessible  in  this  report  makes  the  cyclic  temperature  analysis 
possible  for  composites  with  any  number  of  layers.  Given  the  relevant 
parameters  of  the  heat  flux  intensity,  etc.,  the  computed  temperature 
fluctuations  throughout  the  body  of  the  composite  naturally  lead  to  more 
realistic  determinations  of  internal  thermal  stresses  and  strains. 

3.1  Impulsive  Surface  Heat  Flux 

Even  though  the  development  is  applicable  to  composite  panels 
composed  of  any  number  of  layers,  the  mathematical  development  and 
numerical  presentation  are  limited  to  a  two-layer  configuration  for 
simplicity  and  compactness.  Figure  15  depicts  an  inner  layer,  substrate 
s  or  layer  number  1,  protected  by  an  outer  layer,  coating  c  or  layer 
number  2,  which  acts  as  a  thermal  shield.  With  c  and  s  as  subscripts, 
the  governing  diffusion  equations  in  these  two  regions  are 
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At  the  exposed  surface,  a  constant  heat  flux  Q  is  imposed  impulsively 
on  the  composite  panel,  which  has  an  initial  temperature  of  zero  every¬ 
where.  On  the  back  surface  of  the  substrate  region,  an  insulated 
boundary  condition  is  assumed.  Between  the  two  regions,  the  conditions 
of  equal  temperatures  and  equal  heat  fluxes  are  naturally  valid. 
Altogether,  these  boundary  conditions  are 
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By  using  non-dimensional  variables  defined  by. 


-39- 


=  x/L 


0  »  a  6/L  ' 
s  t 


T  =  k  T  /(OL  )  l 
s  s  s  t  f 


T  =  k  T  /(OL  ) 
c  sc  t 


L  =  L/L 

S  St 


(41) 


the  governing  equations  and  their  associated  boundary  conditions  become 


32t 

3T 

s 

s 

3x2 

96 

a  32T 

9T 

c  c 

c 

a  -s— 2 
s  3x 

96 

3T 

s 

0 

9x 

k  9T 
c  c  _ 

1 

k  nT" 

S  0x 

T  = 

T 

s 

c 

9T  k 

9T 

s  c 

c 

9x  9x 


(x  =  0) 


(x  =  1) 


(X  =  Ls) 


(x  =  L  ) 
s 


(A?) 


(43) 


(44) 


(45) 


(46) 


(47) 


The  group  of  equations  (42)  through  (47)  indicates  that  the  transfent 
solutions  of  Tg  and  Tc  are  dependent  upon  two  parameters:  the  ratio  of 


the  dif fusivities  (a  la  )  end  the  ratio  of  the  conductivities  (k  /k  ) . 

c  s  c  s 
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Finite-Difference  Formulation .  Numerical  solution  of  equations 


(42)  through  (47)  is  most  conveniently  accomplished  by  a  one  time-step 
forward  marching  process  in  conjunction  with  3-point  space  derivatives. 
Although  the  problem  can  be  solved  analytically,  the  resulting  formu¬ 
lation  in  terms  of  segmental  eigenfunctions  becomes  too  unwieldy  to 
handle.  Numerical  treatment  appears  to  be  the  most  expedient  approach. 

As  in  the  case  of  the  spot-heating  problem,  the  ADE-method  (time-splitting) 
is  used  to  insure  speed  and  accuracy.  Based  on  a  uniform  grid  size  of 
Ax,  the  time  step  for  satisfying  numerical  stability  is  governed  by  the 
following : 


A0  <  (Ax)2(as/ac)  (48) 

A0  <  (Ax)  (49) 

which  are  equivalent  to  the  Neuman's  stability  formulation  for  each 
region . 

In  the  computations  carried  out  in  this  analysis,  the  minimum  time 
deduced  from  these  two  criteria  is  further  reduced  by  one  half  to  insure 
accuracy.  The  space  grid  size  Ax  is  determined  by  requiring  the  thinner 
layer  of  the  two  to  have  a  minimum  of  10  divisions.  The  use  of  the 
subscripts  c  and  s  is  dropped,  for  the  different  regions  can  be  distin¬ 
guished  by  referring  to  the  index  of  the  temperature  node  in  the  arrayed 
notations.  By  using  the  ADE-method,  which  is  in  essence  a  split-time 
technique,  equations  (42)  and  (43)  are  replaced  by  a  single  finite-difference 
formula: 
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Figure  15.  Two-Layer  Composite  Configuration  and  Definitions 
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(50) 


(1  +  A) 


where : 


i  =  index  number  of  the  node  (i  =  1,  heating  surface) 

A.  =  (a  /a  )  A  0/  (Ax) 2  ,  in  c-region 
c  s 

A  =  A0/(Ax)2,  in  s-region 

and  subscript  n  denotes  the  value  at  a  later  time-step.  The  node  index 
i  starts  from  the  heating  surface  x  =  1,  i  =  1  and  increases  towards  the 
insulation  surface  of  x  =  0. 

The  heating  boundary  condition  At  x  =  1  (i=l)  is  expressed,  through 
a  threp-point  differentiation  approximation,  by 


T"  =  [4T2  -  T3  +  2(Az)(kg/kc)]/3 


The  interfacial  boundary  condition  is  similarly  obtained  as 


I(4T7  -  Tg)(ks/kc)  +  (4T"  -  t")] 


(52) 


[3  +  3(k  /k  )] 
s  c 


where  the  node  numbers  are  sequential  from  right  to  left  with  node 
number  6  at  the  interface.  The  use  of  these  numerals  is  of  course 
illustrative  only.  Again  the  superscripted  quantities  denote  values  at 
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a  later  time-step,  but  they  become  known  in  the  forward  sweep  algorithm. 
Using  the  same  strategy,  the  end  temperature  at  x  =  0  is  calculated  by: 

t£  =  [4T"  -  t"]/3  (53) 

where  node  number  6  is  used  simply  to  illustrate  the  end  point. 

Results  and  Discussion.  The  calculated  results  are  expressed  in 
terms  of  and  at  the  nodal  points  as  time  proceeds.  Of  these,  the 
most  significant  are  at  the  front,  interface  and  back  positions,  i.e., 
x  =  1,  x  =  L^,  and  x  *  0.  From  a  phenomenological  viewpoint,  all 
temperatures  increase  with  time  and  eventually  become  linearly  dependent 
upon  it.  However,  for  the  purpose  of  analyzing  various  coatings  in 
protecting  the  substrate  layer,  the  relative  magnitudes  of  the  substrate 
temperatures  at  a  fixed  heating  rate  Q  for  a  fixed  duration  0  become  an 
important  criterion  in  coating  selection. 

In  the  numerical  computations  undertaken,  values  of  the  coating 

thermal  conductivity  relative  to  that  of  the  substrate  layer  are  assumed 

as  0.2,  0.5  and  1.  These  are  low-conductivity  coatings  used  as  thermal 

insulators.  For  the  first  set  of  calculations,  the  product  (pC)  of  the 

coating  is  made  equal  to  that  of  the  substrate  layer.  The  relative 

effectiveness  of  the  coatings  can  therefore  be  judged  by  examining  the 

temperature  rises  of  the  substrate  temperatures  at  x  -  and  at  x  -  0. 

The  coating  thickness  is  taken  to  be  0.25  of  the  substrate  thickness. 

Results  of  this  set  of  calculations  are  shown  in  Figure  16,  where  the 

temperature  rises  of  the  substrate  surface  temperatures  at  x  ■  L  are 

s 

grouped  together  for  coatings  of  conductivity  ratios  of  0.2  and  0.5  and 
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1.  The  closeness  of  the  distribution  curves  suggests  that  coatings  with 
low  thermal  conductivities  are  only  mildly  effective  in  protecting  the 
substrate  surface.  The  reason  is,  of  course,  that  in  the  case  of  a 
constant  heat  flux  the  coating  surface  temperature  rises  were  steep  when 
its  conductivity  is  low  and  vice  versa;  the  end  result  is  that  at  the 
interface  with  the  substrate,  the  temperature  changes  with  time  are  not 
materially  different  for  conductivity  ratios  of  0,2  to  1.  The  preceding 
observation  is  not  valid  if  the  heating  environment  is  that  of  con- 
vection. 

Another  set  of  calculations  was  carried  out  with  the  coating 
conductivity  equal  to  that  of  the  substrate  but  with  their  thermal 
capacity  (pC)  ratios  of  1,  1.5,  2,  and  4.  Results  of  varying  the  thermal 
capacity  ratio  indicate  that  the  surface  temperature  of  the  substrate  is 
substantially  reduced  by  using  coatings  with  high  (pC)  values  and  the 
effect  is  much  more  pronounced  than  it  is  when  using  coatings  of  lower 
thermal  conductivities , 

Figure  17  displays  the  front  surface  temperatures  of  the  substrate 
with  coatings  of  various  thermal  capacity  ratios.  The  back  surface 
temperature  variations  with  time  are  shown  in  Figure  18.  Notable 
differences  among  these  response  curves  are  maintained  along  the  time 
scale,  contrary  to  those  in  Figure  16.  Hence  coatings  with  high  rela¬ 
tive  (to  substrate)  values  of  (PC)  are  much  more  effective  than  coatings 
of  low  k-values.  Of  course,  combinations  of  the  two  factors  enhance  the 
effectiveness  of  protection.  Thus  in  relation  to  the  substrate  prop¬ 
erties,  coatings  with  thermal  dif fusivities  lower  than  the  substrates1 
are  preferred  thermal  shields. 
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Figure  17.  Interface  Temperatur  Rises  Due  to  Constant  Heat  Flux  at  Various 
Coating  Thermal  Capacity  Ratios 
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Coating  Thermal  Capacity  Ratios 


3.2  Cyclic  Surface  Temperature 


A  space  structural  component,  when  subjected  to  a  cyclic  surface 
temperature,  may  exhibit  little  internal  temperature  variation  while  the 
internal  temperature  level  follows  almost  in  parallel  with  the  imposed 
surface  temperature  fluctuation.  This  occurs  when  the  surface  tempera¬ 
ture  cycle  is  long  and  the  internal  heat  flow  resistance  is  low.  It 
follows  naturally  that  when  the  surface  variation  is  rapid  in  its 
fluctuation  coupled  with  large  internal  thermal  resistance,  there  would 
be  considerable  internal  temperature  variation  throughout  its  bulk  with 
larger  temperature  excursions  near  the  heating  surface.  In  the  above 
description,  heating  by  cyclic  surface  temperature  variation  or  by 
cyclic  surface  heat  flux  is  considered  identical,  for  one  situation  can 
be  readily  converted  into  the  other  by  the  spatial  gradient  of  the 
former  case. 

The  governing  criterion  for  these  two  extremes  is  the  two  relative 
time  scales:  the  physical  cyclic  time  and  the  thermal  diffusion  time; 
the  former  stems  from  the  imposed  boundary  condition  and  the  latter  is 
an  intrinsic  property  of  the  material  medium.  Compounding  the  situation 
is,  of  course,  the  multi-layer  nature  of  a  composite,  which  makes  the 
analysis  complicated.  From  a  technical  viewpoint,  the  problem  may  be 
stated  as  follows:  given  a  surface  temperature  cycle  history,  determine 
the  steady-state  cyclic  temperature  variation  of  a  multi-layered  compos¬ 
ite  for  which  thermal  properties  are  known  a  priori.  Answers  to  this 


-49- 


question  are  of  significance  in  determining  the  cyclic  thermal  growth  of 
composites  and  sequentially  the  thermal  fatigue  stress  a  component  may 
experience . 

It  turns  out  that  numerical  analysis  by  the  finite-difference  or 
finite-element  method  is  not  the  approach  to  use,  for  the  focus  is  on 
the  steady-state  cyclic  temperature  responses.  By  numerical  methods, 
solutions  would  not  be  forthcoming  until  the  temperature  cycle  becomes 
periodic  and  repeating.  In  principle  it  is  possible;  in  practice  it  is 
not  feasible  to  implement  for  all  kinds  of  combinations  of  problem 
specifications.  Thus  analytical  tools  are  resorted  to. 

Analysis .  Consider  a  multi-layered  composite  depicted  in  Figure 
19.  Starting  from  x  -  0,  the  layers  are  numbered  1  to  n,  each  having 
thermal  properties  distinct  from  other  layers.  At  the  last  layer  (layer 
n) ,  the  exposed  surface  has  a  boundary  condition  described  by  a  cyclic 
temperature  fluctuation  about  a  mean  value.  The  cyclic  part  may  be 
expressed  by 


T  =  £  [a  cosu)  0  +  B  sinoo  0 1  C54) 

s  £=1  L  m  m  m  m  J 

Since  constant  (independent  of  temperature)  thermal  physica]  properties 
are  assumed,  the  entire  problem  is  a  linear  one.  Hence  in  equation  (54), 
only  a  representative  term  needs  to  be  considered;  accordingly,  a 
typical  term  is  detached  from  equation  (54)  and  is  taken  as  the  boundary 
condition  at  x  =  L^.  The  conditions  of  equal  temperatures  and  heat 
fluxes  naturally  apply  at  the  interfacial  positions. 
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Denoting  the  temperature  responses  (from  a  mean  value)  as  ,  T0, 
...  for  these  layers,  the  individual  diffusion  equations  can  be 
expressed  by  an  indexed  notation  as 


9x2 


a.  30 
J 


(55) 


for  j  =  1,  ...  n.  Using  the  total  thickness  L  as  the  reference  length 
and  the  physical  properties  of  the  first  layer  (1)  as  reference,  the 
following  non-dimensional  variables  are  defined: 


x  =  (x/L  ) 
n 

e  =  (a-jO  /Ln2) 

The  governing  equations  for  each  layer  become 

32T.  ax  9T^ 

9x  “j  90 

The  interfacial  positions  are  defined  by 


(56) 


(57) 


x 


v 


L 


n-1 


with  =  1  for  the  exposed  surface,  at  which  the  temperature  of  the 
last  layer  is  specified  by 
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(58) 


^n(x=l)  ^  cos  ^  +  B  sin  uiS 


or 


n (x=l) 


=  (A/2)(eia,e  +  e-ia)0)  -  (iB/2)(eiaJ0  -  e^9) 


(59) 


(Lote  that  the  surface  temperature  variation  in  equation  (58)  has  been 
expressed  in  terms  of  the  non-dimensional  diffusion  time  0  with  a 
corresponding  change  of  cj  to  co  which  will  be  discussed  later.)  To 
determine  the  periodic  solution  for  the  j-th  layer  between  x  =  L  ^  to 
x  =  1^  ,  a  reference  solution  is  obtained  which  satisfies  the  surface 
temperature  fluctuation  (A/2) e^9  .  Let  this  solution  be  called 
A-solution;  then  the  complete  solution  of  the  j-th  layer  due  to  A  coscoQ 
is  the  A-solution  plus  its  complex  conjugate.  To  obtain  the  solution 
due  to  B  sin co©  ,  the  reference  A-solution  is  changed,  (by  replacing  A 
by  -iB)  into  B-solution.  The  complete  solution  due  to  the  surface 
temperature  B  sin  u>0  is  therefore  B-solution  minus  its  complex  conju¬ 
gate. 

A-Solution.  To  obtain  the  solution  T  ,  for  the  j-th  layer,  consider 

•j 

the  substitution 


t  —  /  .  itoG  B  -?  x 

TA__j  -  ( A/2) e  e  3 


(60) 


which  upon  substitution  into  equation  (57)  gives  the  exponential  con¬ 
stant  6.  as 
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Figure  19.  Multi-Layer  Panel  Configuration  and  Definitions 
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Bj  -^(.y a.H 


)  (w/2  V(l+i) 


(61) 


For  the  first  layer  whose  inner  surface  at  x  =  0  is  insulated  against 
heat  conduction,  equation  (60)  can  be  written  as 


TA-1  =  ^A//2^e 


iu)0 


£g^  cosh  |3  x  J 


(62) 


where  is  a  complex  constant.  Note  that  equation  (62)  already  sat¬ 
isfies  the  zero-gradient  requirement  at  x  =  0.  For  layers  from  no.  2 
on,  the  solutions  can  be  expressed  by 

f G.  cosh3.(x  -  L.  )  +  H.  sinh$.(x  -  L#  _  )1 
Lj  J  W  J  ^lJ  (63) 


TA-j  =  (A/2)e] 


j  >  2 


The  functional  form  in  equation  (63)  differs  from  the  more  elementary 
form  of  equation  (60)  for  the  reason  that  equation  (57)  permits  simple 
deduction  of  the  coefficients  G^  and  .  Noting  that  equation  (63) 
applies  for  x  between  L_.  ^  and  L_.  ,  the  starting  value  of  the  terms 
inside  the  brackets  is  simply  G^  which  therefore  equals  the  end  value  of 
the  preceding  layer  no.  (j-1).  In  addition,  since  the  first  term  of 
equation  (63)  has  its  gradient  zero  at  the  starting  position  x  -  L 
the  second  term  has  a  gradient  of  which  must  therefore  be  related 

to  the  gradient  of  the  preceding  solution  at  its  end-point.  In  this  way 
the  calculation  procedure  can  be  made  sequential  and  simple  in  its 
algorithm.  By  proceding  from  the  inner  most  layer  1  to  the  outermost 
layer  n,  the  coefficients  G^  and  H_.  can  be  easily  determined  in  terms 
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of  their  ratios  to  the  first  coefficients  G^.  At  the  surface  position  x 
=  1,  where  the  surface  temperature  (A/2)eia0  must  be  satisfied,  the 
equation  determining  the  first  constant  then  becomes 

G1  _(G^)cosh  V1  -  Ln-1}  +(^)  sinh  en  a  -  Ln  l)J  =  1 


(64) 

where  the  ratios  (G  /G.)  and  (H  /G.)  are  known,  having  been  determined 

n  I  n  I 

by  proceeding  from  layer  to  layer. 

The  A-solution  for  the  (A/2)  eia)0-surf  ace  temperature  perturbation 
is  therefore  embedded  in  equation  (63),  which  contains  complex  coeffi¬ 
cients  G^  and  for  the  j-th  layer  of  the  composite.  Expressing  the 
solution  in  the  form  of 


T  .  =  (A/2)ei,j  [R.  +  il.J 

J  J 


(65) 


where  R_.  and  1^  are  respectively  the  real  and  imaginary  parts  of  the 
complex  function  in  equation  (63) ,  equation  (65)  can  be  expressed  in  a 
more  elementary  form: 

A  r f  f,  r\  T  _  J  ,  .D  \  I  -!  (  T>  _l_  T  n  nr>/.  \P  ^  I  (66) 


T  =  —  Rr  cos  u)0  -  I  sinui0)+  i  (R.sintoB  +  I.  cosoaO) 
A-j  2  L  j  3  J  J 


'Ihe  complete  solution  due  to  a  surface  temperature  fluctuation  of  A 
cos  0)0  is  therefore  the  sum  of  plus  its  complex  conjugate,  and  is 

consequently 


T  =  A  Tr.  cosujG  -  1.  sino)0"l 
A  L  j  j  J 


(67) 
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where  j  is  the  layer  number.  Again,  and  I  are  the  real  and  imagi¬ 
nary  parts  of  the  G-  and  H-terms  in  equation  (63). 

B-Solution.  To  the  surface  temperature  variation  B  sinu)9  ,  there 
is  a  corresponding  temperature  response  for  each  layer.  By  virtue  of 
equation  (59),  the  temperature  response  to  (-iB/2) e^9  is  simply  the 
response  to  (A/2)e  A  if  A  is  replaced  by  (-iB) .  This  is  the 
B-Solution.  Using  equation  (66)  as  the  model,  the  B-Solution  can  be 
expressed  by 


-iB 


B-j 


[<rj 


COSW0  -  I 


j  sinu’e  )+i(Rj  sincoG  +  I  COSCO0)J 


(68) 


Finding  its  complex  conjugate  and  subtracting  it  from  T„  .  gives  the 
temperature  response  function  T^  for  the  j-th  layer  as 


T  =  B[R,  sin  o)0  +  I.  cos  o)0] 
B  J  J 


(69) 


Solution  for  Two-Layer  Composites.  For  a  composite  panel  consist¬ 
ing  of  two  layers,  the  solution  is  not  unduly  complicated  and  can  be 
obtained  by  using  equation  (62)  and  equation  (63)  with  j  -  n  *  2.  At 
the  interfacial  position,  x  =  L^,  T^_^  =  T^  .  Therefore, 


(G2/G1)  =  cosh  61L1 


(70) 


The  equal  heat  flux  condition  at  x  =  L. ,  becomes 
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sinh  3 2L 


(71) 


(h2/g1) 


kl^l 


Introducing  G^G^  and  H^/G^  *nt0  equation  (63)  for  j  *=  2  and  setting 
x  =  1,  an  expression  for  the  fluctuating  outer  surface  temperature  is 
obtained.  The  required  boundary  condition  is  met  by  setting 


G1[(G2/G1)  cosh  B2  (1  -  L1)  +  (H  2/G1)  sinh  &2  ( 1  -  1^)]  =  1  (72) 

The  first  complex  coefficient  (3^  along  with  G 2  and  Ii2>  is  therefore 
determined.  With  the  coefficients  known,  the  temperature  fluctuation 
can  be  separated  into  R  and  I  parts  in  the  region  occupied  by  each 
layer.  The  complete  expressions  for  layers  I  and  2  are  therefore: 


T  =  (A/2)eiUJ0  [cosh  3-xjG. 
A-l  1  J- 


(73) 


1a_2  =  (A/  2)  ela  |^cosh  cosh  32  (x  -  Lp 


— (-  /  7— sinh  3,1*.,  sinh  30  (x  -  L.,  ) 


k2P2C2 


l  l 


']  31 


(74) 


where  G^  is  given  by  equation  (72). 


Dif fuslvity-Cycle  Frequency  co  and  Computational  Results. 


In  the  specification  of  surface  temperature  variation  with  time,  a 
circular  frequency  co  is  used  in  conjunction  with  0,  as  in  equation  (34). 
In  non-dimensional  coordinates,  however,  dif fusi vity-ref erenced  time 
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is  used,  as  defined  by  equation  (56),  which  necessitates  the  use  of  a 


non-dimensional  frequency  co  .  Its  definition  can  be  established  by 


considering  the  identity  co0=co0  ;  hence 


Since  in  the  physical  coordinates,  the  circular  frequency  oa  is  related 
to  the  fundamental  period  P  by  w=  2tt/P,  the  non-dimensional  diffusion-cycle 
frequency  go  is  given  by 


(75) 


which  becomes  a  key  parameter  in  the  analysis  of  periodic  temperature 
responses  of  composites. 

Its  approximate  magnitude  in  a  typical  structural  application  can 
be  established  by  considering  a  low  earth  orbit  with  a  typical  orbiting 
period  of  5A00  seconds  [3].  Used  in  such  a  space  structure  is,  say,  a 
one-inch  thick  graphite-epoxy  composite,  whose  thermal  diffusivity  is 
taken  approximately  as  a  =  3x10  ^  ft2/sec  [4].  With  these  numerical 
values,  u  =  2.7  is  obtained,  thus  establishing  the  range  used  in  this 
analysis . 

First,  for  the  sake  of  demonstrating  the  methodology  and  procedure 
developed,  the  case  of  a  4-layer  composite  is  considered.  All  layers 
are  of  equal  thickness  but  have  thermal  conductivities  in  the  ratios  of 
1,  2,  4,  and  8,  beginning  with  the  inner  layer.  Their  thermal  dif fusivities 
are  taken  to  be  equal.  Figure  20  shows  the  calculated  responses  for  the 
diffusion-cycle  frequency  co  from  1  to  256.  Finite  slope  changes  are 
clearly  indicated  in  Figure  20.  A  general  conclusion  from  Figure  20  is 
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that  at  high  frequencies,  or  at  high  harmonics  of  a  low-frequency 
fluctuation,  thermal  effects  are  confined  to  a  thin  layer  near  the 
surface.  The  frequency-dependent  nature  of  these  distributions  demon¬ 
strates  the  significance  of  the  dif fusivity-cycle  parameter  d>. 

Computed  values  which  made  up  Figure  20  are  found  to  be  most 
conveniently  analyzed  by  the  sequential  method  developed.  If,  instead, 
a  direct  analytical  approach  is  used,  the  necessary  algebra  manipu¬ 
lations  may  become  very  prohibitive. 

Results  for  2-Layer  Composites  (Cyclic  Surface  Temperature) 


Parallel  to  the  computed  results  for  the  case  of  impulsive  surface 
heating  of  two-layered  composites,  representative  calculations  for 
cyclic  surface  temperatures  were  made.  The  outer  layer  (coating)  is 
taken  to  be  0.2  of  the  overall  thickness.  The  first  set  of  calculations 
is  for  reference  only,  with  the  outer  layer  identical  to  the  inner.  For 
the  next  set  of  calculations,  the  outer  layer  is  assigned  a  larger 
thermal  capacity  pC  than  the  inner-layer  value  by  a  factor  in  order  to 
ascertain  the  influence  of  thermal  inertia  of  the  outer  shield.  Subse¬ 
quently,  a  third  series  of  calculations  was  made  with  the  outer  layer 
having  the  same  thermal  inertia  but  with  lower  conductivities  of  0.2  and 
0.1  of  the  substrate  layer.  For  all  the  computations  cited,  the 
dif fusivity-cycle  frequency  co  takes  on  the  values  of  1,  3,  6,  10  and  20, 

thus  covering  an  expected  span  of  variation. 

Figure  21  shows  the  reference  responses  —  in  reality  for  a  single 
layer  composite.  The  in-phase  responses  are  indicated  by  the  R-curves; 
and  the  out-of-phase  responses,  indicated  by  the  I-curves,  are  90°  out 
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of  phase  of  the  cyclic  surface  temperature  variations.  The  I-responses 
therefore  represent  travelling  thermal  waves  or  pulses.  Formation  of 
the  complete  responses  is  supplied  by  equation  (67)  or  (69),  which 
combines  R  and  I  together.  Data  in  Figure  21  suggest  a  demarcation 
criterion  of  d)  =  1.  Below  this  value,  the  internal  temperature  varia¬ 
tions  are  insignificant  and  the  entire  bulk  may  be  treated  as  a  single 
lump  whose  average  temperature  rises  or  falls  in  phase  with  the  imposed 
surface  temperture  cycles.  The  out-of-phase  waves  are  not  insignificant 
to  be  ignored;  for  example,  at  co »  1,  the  back  surface  temperature 
varies  at  0.8  of  the  surface  temperature  excursion  for  the  in-phase 
variation;  however,  the  out-of-phase  variation  has  a  value  of  0.4  of  the 
surface  cyclic  magnitude.  Interpreting  these  temperature  responses  in 
terms  of  the  resulting  thermal  stresses,  the  near  uniform  distributions 
of  the  in-phase  curves  (R)  indicate  that  the  bulk  stresses  rise  or  fall 
with  the  surface  temperature  fluctuations  but  with  only  ±  10  per  vari¬ 
ances  across  the  composite  panels.  Non-uniform  temperature  variations 
and  the  attendant  thermal  stress  variations  lead  naturally  to  thermal 
moments  and  consequently  bending  or  bowing  of  the  panels.  For  the 
out-of-phase  responses,  the  1-curves,  there  is  a  greater  degree  of 
non-uniformity  of  the  temperature  distribution,  ranging  from  zero  at  the 
heating  surface  to  -0.5  at  the  back  surface.  Hence,  the  thermal  bending 
moments  in  the  low-frequency  cases  may  be  due  mainly  to  the  shifted 
temperature  responses . 

As  the  surface  temperature  fluctuation  goes  above  the  critical 
value  of  <£>=  1,  the  alternating  temperature  affect  appears  to  be  more 
confined  in  a  thin  surface  layer,  with  the  remainder  of  the  panel  less 
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affected.  Hence,  for  high-frequency  heating,  damages,  if  any,  are 
confined  to  the  near  surface  region. 

By  taking  the  outer  coat  layer  with  0.2  of  the  overall  thickness 
and  having  a  thermal  capacity  four  times  as  large  as  the  inner  layer, 
the  temperature  responses  were  calculated  and  are  shown  in  Figure  22. 

The  effect  of  larger  (pC)-values  for  the  coat-layer  is  to  have  reservoir¬ 
like  influence  on  the  inner  layer;  the  coat  layer  therefore  alternately 
stores  up  the  in-coming  heat  flux  and  gives  out  what  has  been  stored 
before.  In  cyclic  events,  having  a  coat  layer  with  a  large 
thermal  capacity  is  not  very  influential,  as  compared  with  its  effect  in 
impulsive  heating  analyzed  previously  where  heat  flow  is  not  alternat¬ 
ing.  The  data  in  Figure  22  demonstrate  that  although  there  is  for  both 
in-phase  and  out-of-phase  a  shift  of  the  response  curves  towards  the 
heating  surface,  the  change  is  not  materially  significant  from  the 
curves  in  Figure  21  for  a  single-medium  composite. 

When,  however,  the  outer  layer  has  a  lower  thermal  conductivity 
than  the  inner  layer,  thermal  shielding  thus  afforded  does  substantially 
reduce  the  in-phase  temperature  fluctuation  in  the  inner  layer;  but  less 
so  for  the  out-of-phase  fluctuations.  Data  in  Figures  23  and  24  demon¬ 
strate  respectively  these  phenomena  when  the  outer  layer  has  a  thermal 
conductivity  equal  to  0.2  and  0.1  of  that  of  the  inner  layer.  Particu¬ 
larly  for  high-frequency  heating  u>  =  3  or  above,  the  in-phase  responses 
are  nearly  suppressed  in  the  inner  layer. 

It  is  recalled  that  in  the  case  of  impulsive  heat  load  when  the 
heat  flow  direction  is  not  alternating,  low-conductivity  thermal  shields 
lead  to  very  little  change  in  the  inner  layer’s  temperature-time 
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history.  For  alternating  heat  flow,  however,  this  type  of  shielding 
becomes  particularly  effective. 


3.3  Cyclic  Surface  Heat  Flux 


Having  obtained  composite  panel’s  responses  to  cyclic  surface 
temperature  variations,  the  analogous  problem  with  periodic  surface  heat 
fluxes  is  a  direct  extension  of  the  methodology  established  previously. 
Again,  a  Fourier  series  may  be  used  to  describe  the  surface  heat  flux  of 
which  a  Fourier  component  may  be  represented  by 


q 


q  .  cos  u)0 


qB  sin  u)0 


(76) 


Equation  (76)  is  analogous  to  equation  (54)  for  periodic  surface  temper¬ 
ature  variations.  The  analytical  developments  are  also  parallel  with 
the  previous  case,  up  to  equation  (57).  However,  the  boundary  condition 
at  the  heating  surface  is,  instead  of  equation  (59),  given  by 

k  (3T  /3x)  =  q  A  cos  0)0  4-  qn  sin  o>0  (77) 

n  n  A  & 

Since  the  heat  flux  terms  can  be  expressed  by  the  following  exponential 
terms 

q  =  (q  /  2 )  (e  +  e  )  -  (iqB/2)(e  -  e  )  (78) 


the  complete  solution  is  made  of  individual  solutions,  each  satisfying 
an  individual  heat  flux  component  term  in  equation  (78).  To  obtain  the 
individual  solutions,  a  non-dimensional  temperature  is  defined  by 
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(79) 


TA  “  <klVVA> 


where  is  the  thermal  conductivity  of  the  first  layer,  is  the  tota] 

thickness,  and  q,  is  the  heat  flux  Fourier  coefficient.  For  the  first 
A 

term  of  the  four  in  equation  (78),  the  solution  (A-Solution)  is  assumed 
to  have  the  following  form  for  the  j-th  layer  in  the  composite: 


t  -  / o \  -iu>e  ^jx 

ta_j  -  (l/2)e  e  J 


(80) 


The  exponent,  p,.  is  defined  by  equation  (61).  For  the  first  layer,  the 
solution,  analogous  to  equation  (62),  is  expressed  by 


Vi  ■  (1/2>e 


i  o)6 


cosh 


(81) 


and,  for  other  layers,  by 


( 1/2  )e 


iu)6 


cosh 


3.  (x-L  )  +  H  sinh  3.(x-L  )1 

J  J  1  J  J  J-l J 


(82) 


At  the  heating  surface,  x  =  1,  the  heat  flux  is  to  satisfy  the  first 
term  of  equation  (71);  and  by  using  equation  (75),  the  boundary  condi¬ 
tion  becomes 


SMk/M  [~(G  /G  )  sinh  6  (1-L  )  +  (H  /G  )  cosh  £  (1-L  .  )1  -  1 

niniLnl  n  n-1  n  1  n  n-1  J 


(83) 
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Figure  20. 


Temperature  Functions  of  a  4-Layered  Composite  Due  to 
Periodic  Surface  Temperature  Variation 
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Figure  21.  Temperature  Responses  of  a  Single-Layer  Composite  by 
Periodic  Surface  Temperature  Variation 
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Figure  22.  Temperature  Responses  of  a  Two-Layer  Composite  by 

Periodic  Surface  Temperature  Variation,  (pc)  /(pC)  =  4 
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Figure  23.  Temperature  Responses  of  a  Two-Layer  Composite  by 

Periodic  Surface  Variation,  k  /k  =  0.2 
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Figure  24.  Temperature  Responses  of  a  Two-Layer  Composite  by  Periodic 
Surface  Temperature  Variation,  k^/k  =0.1 
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The  procedure  of  obtaining  the  coefficients  GTs  and  Hfs  is  identical  to 
that  used  previously  for  the  problem  with  a  specified  surface  tempera¬ 


ture  fluctuation.  Equation  (83)  is  used  to  obtain  G^,  instead  of 
equation  (65) . 

With  the  A-solution  thus  obtained  —  defined  by  equations  (81) 
through  (83)  —  which  satisfies  the  heat  flux  boundary  condition  de¬ 
scribed  by  the  first  of  the  four  terms  on  the  right  of  equation  (78), 
the  complete  solution  for  the  surface  heat  flux  q^  cosooQ  can  be 
constructed  by  a  procedure  similar  to  that  for  the  cyclic  surface 
temperature  case,  treated  in  3.2. 

To  be  more  specific,  the  real  and  imaginary  parts  of  the  terms 
inside  the  brackets  of  equation  (82)  are  expressed  by  the  following: 


G.  cosh  $.(x  -  L.  -,)  +  H.  sinh  3 .  (x  -  L.  ') 
J  J  J"1  J  J  3  1 


*  R.  +  il. 

J  J 


(84) 


Then  the  complete  temperature  response  to  the  surface  heat  flux  fluc¬ 
tuation  of  q^  cos  oo9  is  given  by 


-  (k  T  ./L  q  )  =  R.  cos  000  -  I.  sin  oj0 

\-3  1  A-j  n  A  j  j 


(85) 


Similarly,  the  complete  temperature  response  to  the  surface  flux 
q  sinooQ  is 

D 

T„  .  =  (k,T„  ,/L  q  )  =  R.  sin  u)0  +  I .  cos  w0  (86) 

B-j  1  B-j  nMB  j  j 

In  both  equations  (85)  and  (86),  the  layer  index  j  varies  from  1  through 
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n.  An  apparent  structure  of  the  solution  is  that  R_.  is  the  in-phase 
(with  heat  flux)  part  and  1^  the  out-of-phase  part. 


Low  Frequency  Analysis.  Before  discussing  numerical  results  it  is 
important  to  point  out  that  the  out-of-phase  response  1^  is  where  heat 
storage  occurs  and  the  in-phase  response  R_,  is  that  profile  through 
which  heat  conduction  occurs.  To  demonstrate  this  fact,  consider  a 
single-layer  composite.  The  response  to  q  cos  ojg  is  the  following: 

A 

TA-1  =  (klTA-]/LlqA)  =  R1  cos  -  I1  sin  (87) 


where  and  1^  are  derived  from  equation  (81)  and  are  expressed  by 


R^  +  =  cosh  6  x 


j  j^61  sinh 


(88) 


The  complex  coefficient  P1  is  (i  +  l)^ w/2.  The  respective  roles  played  by 
and  I,  in  equation  (87)  can  be  made  clear  by  referring  to  the  dif¬ 
fusion  equation  (57)  for  the  layer.  Integrating  the  diffusion  equation 
from  x*  *  0  to  x  *  1,  where  a  prescribed  heat  flux  of  cos w 9  occurs 
and  using  equation  (87)  for  the  temperature  response,  the  integrated 
form  becomes 


~[ cosuief  R  dx  -  sinu)0  f  T  dxl=  cosa)0 

d0  L  J0  1  J0  1  J 


(89) 


The  right-hand  side  of  equation  (89)  is  the  surface  heat  flux  term;  thus 
equation  (89)  shows  that  heat  flux  is  accumulated  in  the  I^-term  and 
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correspondingly  the  boundary  condition  of  heat  flux  at  x  =  1  is  sat¬ 
isfied  via  the  R^-term. 

Because  of  the  above  examination  of  the  terms1  composition,  it  is 
of  importance  to  note  that  if  the  cyclic  heat  flux  condition  is  such  that 
heating  or  cooling  takes  place  over  a  long  period  cf  time,  i.e.,  cj  -*  0 
or  (3  ->0,  then  the  heat  accumulation  terms  I's  tend  to  be  large  since 
there  is  more  time  to  pile  up.  To  substantiate  that,  the  one-term 
solution  of  equation  (88)  can  be  decomposed  into  its  real  and  imaginary 
parts:  Since  (3  is  a  complex  constant,  asymptotically  small  values  of  (3 

or  cj  lead  to  the  following  series  expansion  in  terms  of  cj  as  the 
parameter : 


cosh  Bx 

8  sinhS 


(90) 


Thus  the  imaginary  part  T,  is  inversely  proportional  to  cj  or  directly  to 
the  cyclic  period,  a  clear  indication  of  the  heat  accumulation  effect. 
Moreover,  for  quasi-steady  state  heating  or  cooling,  i.e.  ,  co  0  the 
real  part  R,  becomes  asymptotically  parabolic,  as  indicated  by  the 
leading  term  of  equation  (90) . 

Tt  is  clear  that  for  any  combination  of  the  physical  parameters, 
whether  it  is  a  single-layer  or  multi-layer  composite,  there  invariably 
exist  asymptotic  distributions  of  the  in-phase  and  out-of-phase  parts,  R 
and  I  for  different  layers.  Their  mathematical  representations 
may  be  more  involved,  but  the  existence  can  be  easily  established. 

Numerical  Results  for  2-Layer  Composites  (Cyclic  Flux).  As  a 
reference  combination,  consideration  is  given  to  a  two-layer  composite 
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with  both  layers  having  identical  physical  properties.  Tn  essence,  it 
is  a  single  layer  analysis.  Computed  response  functions  R  and  I  are 
presented  in  Figure  25  for  diffusion-cycle  frequencies  of  <I>  =  1,  6,  10 
and  20. 

It  is  clear  from  the  graphical  display  that  for  0=  1,  the  dis¬ 
tribution  of  I  is  almost  asymptotic,  which  is  given  by  equation  (90)  as 
-1.  A  slight  variation  from  x  «  0  to  x  *  1  is,  however,  noted,  indicat¬ 
ing  its  dependency  on  go  yet.  For  the  real  part,  the  distribution  for  co 
=  1  is  almost  indistinguishable  from  that  in  equation  (90)  for  co  ^  0. 
Translating  these  observations  into  thermal  stress  considerations,  it  is 
the  bulk  temperature  rise  that  yields  thermal  stresses.  Hence  it  is  the 
I-distribution  that  governs  the  stress  magnitudes;  and  at  low  frequen¬ 
cies  there  would  be  higher  stresses.  The  thermal  stresses  caused  by  the 
R-distributions  are,  however,  in  self-equilibrium  and  these  stresses 
contribute  a  thermal  moment  across  the  composite  slab,  resulting  in  its 
flexural  bending  as  a  consequence.  Naturally,  as  the  frequency  increas¬ 
es,  the  temperature  distribution  functions  R  and  I  diminish  In  magnitude 
and  so  do  the  thermal  stresses.  These  qualitative  trends  also  apply  to 
multi-layer  composites. 

With  two-layer  composites,  numerical  computations  were  completed 
first  with  the  outer  layer  (0.25  of  the  inner  layer  thickness)  having  a 
larger  thermal  capacity  of  4  times  the  inner  value.  The  external 
coating  behaves  as  a  thermal  sink  with  respect  to  the  surface  heating 
and  results  in  a  general  reduction  of  the  temperature  responses  in  the 
inner  layer.  Figure  26  indicates  the  response  curves  for  R  and  I  in  the 
same  frequency  range  as  in  Figure  25.  While  the  real  part  R  (in-phase 
with  surface  flux)  shows  very  little  change  from  those  in  Figure  25,  the 
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imaginary  part  I  (out-of-phase  and  for  heat  storage)  indicates  signifi¬ 
cant  reductions  in  magnitude,  thus  reflecting  the  influence  of  the 
thermal  capacity  in  the  outer  layer. 

The  other  parameter  of  interest  is  the  thermal  conductivity.  By 
decreasing  the  coat  layer  conductivity  to  0.2  of  the  inner  layer  (with 
equal  therma]  capacity  for  both),  the  response  curves  are  shown  in 
Figure  27.  A  notable  feature  is  that  the  1-distributions  are  very  much 
similar  to  those  in  Figure  25,  indicating  little  effect  of  low  conducti¬ 
vities  of  the  coat  layer  on  the  level  of  the  bulk  temperature  change. 
There  are,  however,  considerable  increases  of  the  R-distribut ions  over 
those  of  Figure  25  for  single  layer  performances.  With  the  outer  layer 
acting  as  a  heat  barrier,  temperature  differentials  across  the  barrier 
are  increased  to  account  for  the  specified  surface  heat  flux.  Hence 
from  a  thermal  protection  viewpoint,  low  conductivity  in  the  outer  layer 
terds  to  produce  greater  thermal  stresses. 


4.  CONCLUSIONS 


In  this  report,  two  transient  heat  conduction  problems  for  applica¬ 
tion  to  the  thermal  evaluation  of  composite  materials  are  analyzed  and 
discussed . 

The  purposes  of  these  two  analyses  are  at  least  two-fold:  One  is 
to  demonstrate  by  detailed  numerical  results  that  composite  materials 
with  a  larger  in-plane  thermal  conductivity  than  the  transverse  value 
exhibit  lower  surface  temperatures  when  the  thermal  loading  consists  of 
irradiation  by  a  concentrated  cylindrical  beam.  Fffective  radial  heat 
spread  along  the  in-plane  direction  serves  to  reduce  local  heat 
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accumulation  and  thereby  mitigates  heat-spot  temperature  rises.  Accord¬ 
ingly,  preventing  or  delaying  damage  due  to  high-intensity  thermal 
radiation  on  composite  surfaces  can  be  enhanced  by  using  materials  with 
large  in-plane  thermal  conductivities.  The  second  objective  of  the 
investigation  is  to  develop  a  methodology  whereby  composite  panels 
undergoing  cyclic  heating  and  cooling  can  be  analyzed  for  their  periodic 
temperature  responses.  This  is  particularly  important  for  determining 
thermal  stresses  resulting  from  alternating  temperature  fluctuations. 
The  methodology  developed  is  especially  use.ful  since  numerical  approach¬ 
es  can  be  very  time-consuming  and  unreliable. 
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Figure  25.  Temperature  Responses  of  a  Single-Layer  Composite  by 
Periodic  Surface  Heat  Flux 
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Figure  26.  Influence  of  Large  Coating  Layer  Thermal  Capacity  on 
Temperature  Responses  of  a  Two-Layer  Composite  by 
Periodic*  Surface  Heat  Flux 
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Figure  27.  Influence  of  Low  Coating  Layer  Thermal,  Conductivity  on 
Temperature  Responses  of  a  Two-Layer  Composite  by 
Periodic  Surface  Heat  Flux 


-7  7- 


KFFhKFNCES 

].  Roaohe,  P.J.,  Computational  Fluid  Dynamics,  llcrmosa  Publisher:-. 
1976,  pp.  95-%. 

7.  Cars  lav  9  H.S.,  and  Jaeger,  J.C.,  Conduct  1  on  r£  lira  t  in  Solid.'  , 
Oxford:  Clarendon  Press,  1947. 

3.  Thorton,  F.A.,  and  Paul,  D.B.  ,  nThe rma 1-St rue  turn!  Analysis  ui 
large  Space  St rue tures , 11  ATAA  Paper  No.  63-1018,  to  be  published 
in  AIAA  Journal . 

Fan,  L.  and  Eoyee,  W.F.,  "Thermal  Conductivities  and  Di f fusivi t ie*. 
of  Graphite- Fpnxy  Composites,"  AFWAL-TR-83-3002 ,  Feb.,  1983. 


-78- 


NOMENCLATURE 


a 

C 


i.j 


k 

L 

P 


Q»q 


Vq 


B 


r 

T 

w 

X 

z 

R 

I 


heat  beam  radius 

specific  heat 

node  indexes 

thermal  conductivity 

layer  thickness 

period  of  fluctuation 

heat  flux  intensity 

Fourier  coefficients  for  heat  flux 

radial  coordinate  from  spot  (beam)  axis 

temperature  rise 

panel  thickness 

depth  coordinate  from  insulated  surface 
depth  coordinate  from  heating  surface 
real  part  of  a  complex  function 
imaginary  part  of  a  complex  function 


Mon-dimensional  Variables 


z  =  (z/a) 

r  =  (r /a)  k  /k 
r  z 

Q  =  (a  6/a2)  or  (a  0/L  2) 

7  St 

T  =  (k  T/aQ) 
z 
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Greek  Symbols 


A 

6 

a 

P 

w 

P 


Subscripts 


b 

c 

j 

o 

r 

s 

t 

z 


grid  size  A  =  Ar  =  A  x  =  A  z 
nondimensional  time  step 
thermal  diffusivity  k/pC 

density 

circular  frequency 
frequency  parameter 


back  surface,  beam  center 

coating  layer 

ordinal  number  for  layer 

beam  spot  center 

radial  direction 

substrate  layer 

total 

depth  direction 
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APPENDIX 


Comparison  of  Results  with  Different  Grid  Sizes 

It  is  well  established  that  computational  results  based  on  the 
finite-difference  approach  are  sensitive  to  the  grid  size  used.  The 
finer  the  grid  size,  the  more  accurate  are  the  results;  but  far  more 
time-consuming  is  the  task.  Halving  the  grid  size  usually  increases  the 
process  time  by  a  factor  of  eight  and  more.  For  this  reason,  some 
balance  is  needed  between  accuracy  and  excessive  computational  effort. 

In  the  computational  effort  undertaken  is  this  report,  the  grid 

size  is  obtained  by  dividing  the  non-dimensional  heat  spot  radius  into 

five  divisions.  Since  the  non-dimensional  radius  is  defined  by 

r  =  (r/a)Jk  /k  where  a  is  the  beam  radius,  each  division  of  r  is 
*  z  r 

therefore  equal  to  (ykT/k^)/5.  In  the  depth  direction,  the  dimension¬ 
less  variable  is  z  =  (z/a),  each  step Az  is  made  equal  to  Ar,  and  the 
number  of  divisions  is  determined  by  the  width  w  and  the  ratio  of  k  /k  . 
Tn  this  way,  a  grid  uniform  in  the  non-dimensional  coordinates,  but 
non-uniform  in  the  physical  coordinates,  is  achieved. 

In  order  to  ascertain  typical  differences  in  the  results  due  to 
different  grid  sizes,  the  case  of  k^/k^  *  1  and  w/a  =  1  was  analyzed  by 
two  parallel  computations:  one  using  5  divisions  and  the  other,  10 
divisions  for  the  heat  spot  radius.  Results  of  these  two  comparative 
calculations  are  presented  in  Figure  A-l,  which  shows  the  spot  center 
temperature  rises  as  heating  proceeds.  Naturally,  the  results  based  on 
10  divisions  for  the  heat  spot  are  more  reliable  than  those  based  on  a 
coarser  grid.  However,  the  difference  in  the  temperatures  at  the  end  of 
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the  non-dimensional  time  0-  (c^0/a2)  =  A  is  quite  small  compared  to 
their  mean  value:  the  difference  being  0.05  out  of  theii  average  value 
of  1.2  at  the  spot  center.  At  the  back  surface  of  the  slab,  opposite  to 
the  spot  center,  the  difference  becomes  0.07  out  of  their  mean  value  of 
0.9.  The  computer  process  times  were  A  seconds  and  55  seconds  respecti¬ 
vely.  For  other  conductivity  ratios  k^/k^<  1  with  a  much  finer  grid  in 
the  z-direction  than  in  the  r-direction,  the  computer  process  time  was 
found  to  be  much  more  than  55  seconds. 

Additional  comparisons  of  these  two  parallel  sets  of  results  are 
presented  in  Figures  A- 2  and  A-3.  The  former  shows,  at  the  end  of  heating 
time  0  =  10,  the  isothermal  lines  of  temperature  0.8,  0.6  and  0.A  of  the 
spot  center  temperature.  The  results  using  different  grid  sizes  are 
nearly  coincidental  to  each  other.  Thus,  the  use  of  a  smaller  number  ol 
division  is  deemed  adequate.  Shown  in  Figure  A-3are  the  temperature 
profiles  in  the  r-direction  on  the  front  and  back  surface  at  0  =  10. 

Even  though  the  spot-center  temperature  rises  are  1.A9  and  1 . A 1  for  5 
and  10-direction  calculations,  their  normalized  (with  respect  to  the 
spot-center  temperature)  curves  are  parallel  to  each  other.  Thus,  the 
use  of  5-division  appears  quite  adequate  for  the  demonstrative  analysis 
presented  in  this  report. 
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Figure  A-l.  Effect  of  Grid  Size  on  Temperature  Calculations,  Spot  Hearing 
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Figure  A-2.  Effect  of  Grid  Size  on  Isothermal  Contours,  Spot  Heating 
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Figure  A-3.  Effect  of  Grid  Size  on  Surface  Temperature  Profiles, 
Spot  Heating 
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